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ABSTRACT 


Grand Canonical Monte Carlo simulations are used to reproduce the N 2 /CO 
ratio ranging between 1.7 x 10 -3 and 1.6 x 10 -2 observed in situ in the 
Jupiter family comet 67P/Churyumov-Gerasimenko by the ROSINA mass 
spectrometer aboard the Rosetta spacecraft, assuming that this body has been 
agglomerated from clathrates in the protosolar nebula. Simulations are done 
using an elaborated interatomic potentials for investigating the temperature 
dependence of the trapping within a multiple guest clathrate formed from a 
gas mixture of CO and N 2 in proportions corresponding to those expected 
for the protosolar nebula. By assuming that 67P/Churyumov-Gerasimenko 
agglomerated from clathrates, our calculations suggest the cometary grains must 
have been formed at temperatures ranging between ~31.8 and 69.9 K in the 
protosolar nebula to match the N 2 /CO ratio measured by the ROSINA mass 
spectrometer. The presence of clathrates in Jupiter family comets could then 
explain the potential N 2 depletion (factor up to ~87 compared to the protosolar 
value) measured in 67P/Churyumov-Gerasimenko. 


Subject headings: astrobiology - comets: general - comets: individual 
(67P/Churyumov-Gerasimenko) - solid state: volatile - methods: numerical 
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Introduction 


The thermodynamic conditions prevailing in many bodies of the solar system suggest 


that clathrates could exist in the Martian permafrost (Thomas et al. 2009; Swindle et al 


2009 Mousis et al.||2013 ), on Titan ( Choukroun fe Sotin||2012 Mousis et al.||2014 ), as well as 
in the interiors of other icy satellites (Kieffer et al. |2006 ; Hand et al. 2006). It has also been 
suggested that the activity observed in some cometary nuclei results from the dissociation 


of these crystalline structures (Marconi & Mendis 1983; Smoluchowski 1988; Marboeuf et 


al. 2010, 2011, 2012 Mousis et al. 2012). Several indirect evidences suggest that clathrates 


probably participated in the formation of planetesimals and the building blocks of giant 


planets in the outer solar system (Lunine & Stevenson 1985 Mousis et al. 2010, 2014) 


In the absence of experimental data existing at low-temperature (20-200 K) and 
low-pressure conditions (10~ 13 -10“ 3 bar), which are typical of those encountered in 


planetary environments (Lunine & Stevenson 1985 Sloan & Koh 2008), clathrates are 


characterized via theoretical modeling. The approach usually employed in planetary science 


is based on the statistical mechanics model initially proposed by van der Waals & Platteeuw 


(1959), which makes use of simplified intermolecular potentials calibrated on equilibrium 


measurements performed at relatively high temperatures. In consequence, the use of these 
potentials at thermodynamic conditions relevant to those of the protosolar nebula (hereafter 
PSN) for predicting the composition of clathrates deserves to be confronted with more 
sophisticated approaches. 

In this present work, we aim at reproducing the N 2 /CO ratio ranging between 1.7 
x 10~ 3 and 1.6 x 10~ 2 observed in situ in the Jupiter family comet 67P/Churyumov- 
Gerasimenko (hereafter 67P) by the ROSINA mass spectrometer aboard the Rosetta 


spacecraft (Balsiger et al. 2007), which is found to be depleted by a factor up to ~87 


compared to the value of 0.148 hypothesized for the protosolar nebula (see below). By 
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assuming that 67P has been agglomerated from clathrates, it is possible to derive the 
temperature range of formation of these crystalline structures in the PSN by mean of Grand 
Canonical Monte Carlo (GCMC) simulations based on elaborated interatomic potentials. 
These allowed us to investigate the temperature dependence of the trapping within a 
multiple guest (hereafter MG) clathrate formed from a gaseous mixture of CO and N 2 in 
proportions corresponding to those expected for the protosolar nebula. 


2. Computational details and simulation procedure 


We assume multiple guest (MG) clathrate formation from a gaseous mixture composed 
of N 2 and CO in proportions specified from a plausible protosolar gas phase composition. 
For this, we have assumed that both all C and N are in form of CO and N 2 in the initial 
disk’s gas phase. This assumption is in agreement with thermochemical models of the 


protosolar nebula (Lewis & Prinn 1980 Prinn & Fegley 1989 Mousis et al. 2002). The use 


of C and N protosolar elemental abundances from the compilation of Lodders et al. (2009) 


allowed us to derive a gaseous mixture with mole fractions of 0.871 for CO and 0.129 for 
N 2 , giving a N 2 /CO ratio of 0.148 in the PSN. 

The MG clathrate composition has been computed along its equilibrium pressure curve 
P^,mg given by (Thomas et al. 2007 Mousis &; Schmitt| 2008): 


P 


eq,MG 


-1 


E^ 

. 1 eq,i 


( 1 ) 


where iji is the mole fraction of species i in the gas phase (here CO or N 2 ) and P eq ,i the 
equilibrium pressure of single guest clathrate formed from component i only. 


The equilibrium pressures P cqi j (in bar) are determined by using an equation based on 


an Arrhenius law (Miller 


1961): 
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In P eq!i = ^ + Bi, 


( 2 ) 


where A* and Bi are constant parameters depending on the nature of the species trapped 
in the clathrate and T is the temperature (K). A j and B f have been determined by fitting 


the available theoretical and laboratory data (Mousis et ah 2008) and their values are given 


in Table [I] The calculated values of P eqi ; and partial pressures of CO (Pco) and N 2 (Pn 2 ) 
are represented as a function of the inverse temperature in Fig. [T[ Here, the equilibrium 
pressure of MG clathrate is in the ~5.2 xlO~ 10 -2.9 xl0~ 3 bar range when T is varied 
between ~52 and 100 K, respectively. 

We have calculated the composition of N 2 -CO clathrates via Monte-carlo (MC) 


simulations in the Grand Canonical ensemble (GCMC) (Frenkel & Srnit 2002) for 


temperatures ranging from 52 to 100 K (the choice of this temperature range is explained 
in Sec. [3j but it is worth noting that, below 50 K, the equilibration time is in fact too long 
for the GCMC simulation), with an increment of 4 K between each computation. All our 
calculations have been performed in the case of Structure I (SI) clathrates (see structural 


details in Sloan & Koh (2008)). This choose has been motivated by the fact that CO is the 


dominating species in the gaseous mixture and is known to form Structure I single guest 


clathrate (Mohammadi 2005). 


In our system, the considered crystal size consists in 125 cubic unit cells (5x5x5), 
corresponding to 5,750 water molecules. The dimension of one parameter of the cubic 
simulation box is set equal to 60.15 A during all our simulations. Periodic boundary 
conditions are applied to mimic infinite crystal. The water molecules are modeled using 


the well-known TIP4P/2005 model (Abascal & Vega 2005), allowing them to translate and 


rotate during the simulation. Models for N 2 and CO molecules are taken from Potoff & 


Siepmann (2001) and Piper et al. (1984), respectively. One hundred million MC steps were 
























6 


performed including insertion, deletion, translation and rotation of the molecules. Only the 
last 50 million steps were used to compute the data. The first 50 million steps have been 
discarded from the analysis and were only used to equilibrate the system. 

From the knowledge of the partial pressure of each species (see Fig. [Tj), GCMC 
simulations have been performed to compute the composition of the MG clathrate. 
However, the number of MC steps needed to equilibrate the system strongly increases 
when temperature decreases. Simulation tests showed that below a temperature of 52 K, 
this number becomes even larger than the 50 million steps we simulated. At these low 
temperatures, statistically relevant results on the mole fractions of encaged molecules in 
our MC calculations thus become too time consuming. In consequence, below 52 K, these 
quantities have been evaluated via a thermodynamic extrapolation, described below, and 
fitted to the results computed at higher temperatures. 

After the equilibration period, chemical equilibrium takes place between gas and 
clathrate. It is then possible to write an equilibrium constant K t for each species i in the 
form: 


Ki = 


Pi/P° 


( 3 ) 


where a* is the activity of species i in clathrate defined from the mole fraction of species 
i in clathrate (xi) and its activity coefficient (y; ) as a* = y^ay. The evolution of K t with 
temperature should obey the Van’t Hoff relation: 


dln/y AEi 
d T ~ RT2’ 

where A E{ is the entrapping energy of CO or N 2 , and R the ideal gas constant. 


Under the simulated conditions, the occupancy of the clathrate cages is still maximum. 





7 


We thus, obtained an amount of 8 molecules of CO and N 2 per unit cell, implying 
that the total number of molecules trapped in clathrate is constant, independently of 
the temperature. By Neglecting the non ideal terms of the activity coefficient at first 
approximation, this allowed us to get 7 * = 1 and to rewrite Eq. [| as follows: 


din (TV-/p.) _ A Ei 

d(l/T) R ’ [ } 

with N t the number of encaged molecules of type i. Our simulations allowed retrieving the 
values of Nqo and TVn 2 at temperatures higher than 52 K and provided us with a direct 
access to the different values of the quantity In(TVj/Pj). Figure [ 2 ] represents these quantities 
plotted for CO and N 2 as a function of 1 /T. They show a linear behavior with a correlation 
coefficient higher than 0.999 for each species, in agreement with Eq. [5] The computed data 
have been fitted via the use of linear equations, allowing to find ln(iVco/-Pco) = 1685.145/T - 
15.617 and ln(iV N 2 /P N2 ) = 1554.469/T - 15.972. The entrapping energies comes directly 
from the fit, APco = — 1685.145xR, and A — — 1554.469xR. These two linear equations 
have been used to estimate the Nm/Nco ratio at temperatures lower than 52 K. Note that, 
for simplification, we will use below the abbreviation N 2 /CO for this ratio. 


3. Results 

Figure |3] shows the evolution of N 2 /CO ratio as a function of temperature in the 
20-100 K range. This ratio monotonically increases with the growing temperature. The 
figure exhibits a linear regime in the 50-100 K range whereas at lower temperatures, the 
curve trend is less steep and the ratio converges smoothly towards zero. The N 2 /CO ratio 
found is equal to ~1.5 x 10 ~ 4 at 20 K. This value is 50 times smaller than at 50 K. In the 
temperature range considered here, the calculated N 2 /CO ratio is significantly lower than 




the ratio of ~0.15 in the coexisting gas phase. This indicates that the clathrate formation 
favors the CO entrapping at the expense of N 2 . This behavior is related to the difference in 
entrapping energy between CO and N 2 . Indeed, AAco being lower, the entrapping of CO 
is selectively favored when the temperature decreases. 


Figure [4] is a zoom of a portion of Fig. [3] given for easy reading of the correspondence 
between the N 2 /CO ratio measured in 67P by the ROSINA instrument and the formation 
temperature of the ice grains from which the comet agglomerated. Taking into account the 
strong variation of the N 2 /CO measurement between 0.17 to 1.6% depending on the position 
of the Rosetta spacecraft above the surface of the comet nucleus (Rubin et al. 2015), we find 
that the ice grains at the origin of 67P formed at temperatures ranging between ~31.8 and 
69.9 K in the protosolar nebula, with corresponding equilibrium pressures ranging between 
6.0 x 10~ 19 and 2.1 x 10 -6 bar. For the sake of information, the mean N 2 /CO ratio of 


0.57% corresponding to the averaging of the 138 spectra obtained by Rubin et al. (2015) is 
represented on Fig. [4j The corresponding formation temperature of the ice grains is of ~45 
K in the PSN, with an equilibrium pressure of 3.3 x 10~ 12 bar. 


4. Discussion and Conclusions 

The composition of a MG clathrate formed from a gaseous mixture of N 2 and CO 
in proportions corresponding to those expected for the protosolar nebula (87.1 % for CO 
and 12.9% for N 2 ) has been investigated in the 20-100 K temperature range. Above 
50 K, the clathrate composition has been computed via Grand Canonical Monte-Carlo 
simulations for pressures ranging from ~5.2 xlO -10 to 2.9 xlO -3 bar. Below 50 K, the 
clathrate composition has been extrapolated via the use of a Van’t Hoff relation. The 
results show that, at thermodynamic conditions relevant to those of the protosolar nebula, 
CO has a much higher propensity than N 2 to be trapped in clathrates. Assuming that 
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67P agglomerated from clathrates, our calculations suggest that the cometary grains must 
have formed at temperatures ranging between ~31.8 and 69.9 K in the protosolar nebula to 
match the N 2 /CO ratio measured by the ROSINA mass spectrometer (Rubin et al. 2015). 

Whilst narrower, the range of formation temperatures inferred from our model for the 
grains of 67P is consistent with the one (~22-80 K) found from the reading of Fig. 2 of 


Mousis et al. (2012) who performed calculations of planetesimals compositions based on the 


classical statistical mechanics model of van der Waals & Platteeuw (1959). In the absence 
of experiments at this temperature range, the fact that these two different approaches lead 
to similar conclusions, namely that clathrates can explain the N 2 /CO ratio observed in 67P, 
suggest that this scenario is plausible. The presence of clathrates in Jupiter family comets 
could then explain the apparent N 2 depletion (factor up to ~87 compared to the protosolar 
value) measured in 67P. 
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Table 1: Parameters of the equilibrium curves of the considered single guest clathrates 


Molecule type i 

A / (K) 

Bi 

CO 

-1685.54 

10.9946 

n 2 

-1677.62 

11.1919 



Pressure (Pa) 
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Fig. 1.— Calculated pressure and partial pressures of CO and N 2 in the gas as a function of 
the inverse temperature based on the Arrhenius law and using parameters in Table [l] 
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Fig. 2.— Evolution of the number of encaged molecules in the cases of CO (a) and N 2 (b) as 
a function of the inverse temperature. A linear fit was performed on the data with correlation 
coefficients higher than 0.999 in both cases (see text). The figure shows the results of the 
GCMC simulation between 52K and 100K. 


i,gas 
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Fig. 3.— N 2 /CO ratio in clathrate as a function of formation temperature. The results are 
derived from the GCMC simulations performed above 52K. Below this temperature, they 
are based on the entrapping energies derived from the GCMC simulations. 
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Fig. 4.— Minimum and maximum N 2 /CO ratios measured in 67P and corresponding for¬ 
mation temperatures for the ice grains. The results are derived from the GCMC simulations 
performed above 52K. Below this temperature, they are based on the entrapping energies 
derived from the GCMC simulations. 







